Numerical simulations of a two-dimensional 
lattice grain boundary model 



A. Jaster, H.H. Hahn 

Institut fiir Theoretische Physik, TU Braunschweig, 
Mendelssohnstr. 3, D-38106 Braunschweig, Germany 

May, 1997 



Abstract 

We present detailed Monte Carlo results for a two-dimensional grain boundary 
model on a lattice. The effective Hamiltonian of the system results from the micro- 
scopic interaction of grains with orientations described by spins of unit length, and 
leads to a nearest-neighbour interaction proportional to the absolute value of the 
angle between the grains. Our analysis of the correlation length £ and susceptibility 
X in the high-temperature phase favour a Kosterlitz-Thouless-like (KT) singularity 
over a second-order phase transition. Unconstrained KT fits of x an d £ confirm 
the predicted value for the critical exponent u, while the values of r\ deviate from 
the theoretical prediction. Additionally we apply finite-size scaling theory and in- 
vestigate the question of multiplicative logarithmic corrections to a KT transition. 
As for the critical exponents our results are similar to data obtained from the XY 
model, so that both models probably lie in the same universality class. 



1 Introduction 



Phase transitions in two-dimensional systems with continuous symmetry are of great 
interest since many years. In these systems conventional long-range order at non-zero 
temperature cannot occur JEf. However, two-dimensional particle systems or the XY 
model are characterized at low temperatures by quasi-long-range order, while the high- 
temperature phase is disordered. The critical behaviour of the models depends on the 
dimension of the order parameter. For the two systems mentioned the order parameter 
dimension is two. Physical realizations of two-dimensional 0(2) symmetric systems are 
films of liquid helium or superconductive layers. 

The XY model provides a simple model of a system with continuous symmetry. For the 
phase transition numerical studies favour a Kosterlitz-Thouless-like [|] behaviour. The 
KT mechanism is based on a topological defect (called vortex) unbinding scenario. In the 
low-temperature region vortices are bound in pairs. The dominant excitations are spin 
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waves, which destroy longe-range order. Increasing the temperature leads to an unbinding 
of vortex pairs. The theory predicts a continuous transition from the quasi-long-range 
phase to the disordered phase derived from renormalization group treatment. 

Particle systems are characterized by two order parameters related to different topo- 
logical defects. Therefore, for the melting transition additional complications arise. There 
are several theoretical approaches for this transition. Halperin and Nelson [Q] as well as 
Young || enhanced the ideas of Kosterlitz and Thouless. The KTHNY theory deals with 
a dislocation unbinding mechanism, which is responsible for the melting transition, and 
a disclination transition, which destroys the nearest-neighbour-bond orientation. An al- 
ternative scenario has been proposed by Chui ||. He presented a theory via spontaneous 
generation of grain boundaries, i.e. collective excitations of dislocations. He found that 
grain boundaries may be generated before the dislocations unbind if the core energy of 
dislocations is sufficiently small, and predicted a first-order transition. 

Another possibility to study two-dimensional melting is the simulation of the defect 
system itself. Saito performed Monte Carlo simulations of dislocation vector systems |5| 
and found that for systems with large dislocation core energy unbinding of dislocation 
pairs causes a continuous phase transition, while small core energies produce a first-order 
transition by formation of grain boundaries. 

Patashinskii f?J proposed a model, which can be seen as a mesoscopic lattice model 
for grain boundaries. The energy per unit length of a boundary between two grains 
with small angle A0 = — is proportional to the absolute value of A<f>. A simple 
ansatz for the description of the interaction of grains is then given by putting them on 
a square lattice with nearest-neighbour interaction proportional the absolute value for 
small lattice gradients. The orientation of the grains is described by spins of unit length 
s(x) = (cos0(x), sin0(x)). Therefore Patashinskii used a Hamiltonian of the form 



The number N is related to the symmetry of the system. For example, particle systems 
with a hexagonal symmetry are described by N = 6. However, the partition function is 
essentially independent of N. Therefore we neglected this parameter and used angles of 
range < 0(x) < 2tc. Also, we used the simplified model 



for numerical simulations, since (at least in the case of a continuous phase transition) 
only the behaviour for small A0's is relevant. This coincides with an XY model with the 
negative cosine replaced by its absolute value. In this paper we try to determine the order 
and mechanism of the phase transition of this model. 

The KT scenario predicts an exponential singularity for the correlation length 
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if t — /3 C — P \ 0, while the specific heat Cy should not show any divergent behaviour. 
The critical exponent r\ defined by 



x ~ e~\ (5) 

and the critical exponent v are given by 

V = 2 ' V = 4 ' 

while is a non- universal constant and b x = (2 — 77)6^. An alternative approach is a 
conventional second-order behaviour with power-law singularities of 

C(t) = a^t~ v (7) 

and 



X{t) = a x t 
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In the following we analyze our extensive Monte Carlo results to answer the question of 
kind and order of the phase transition of the grain boundary model. 



2 Simulations 

2.1 Algorithms and measurement 

For the simulations we used a combination of the Metropolis, over-relaxation M, 01 and 



cluster |10| , algorithm, depending on the temperature and size of the lattice. The 
study was carried out on lattice sizes of L 2 = 60 2 , 120 2 , 240 2 , 480 2 and 960 2 with periodic 
boundary conditions. 

For each simulation we measured the energy, specific heat, 'magnetization', susceptibil- 
ity, fourth-order cumulant and the zero momentum correlation function (for a definition 
see eq. (0) below). In doing so, we used conventional estimators for all observables. 
Additional improved estimators |T2[ were measured for the susceptibility and correlation 



function. Typically, errors were of the order of 1% or less. 

For the calculation of the specific heat we used two different ways: on one hand we 
measured the fluctuations in the energy per point 

C v = p 2 L 2 ((e 2 )-(e) 2 ) (9) 

on the other hand we compute the derivative 

C y = -(3 2 de/d(3 . (10) 

We find that both methods give consistent results as shown in fig. |l|. Except for the 
smallest lattice the data do not show any significant finite-size effects. As one can see 
the specific heat develops a smooth peak at an inverse temperature below the transition 
point (3 C ~ 0.83. At the transition point the energy and the specific heat stay finite and 
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Figure 1: Specific heat as a function of (3. Full symbols stand for the calcu- 
lation by numerical differentiation, open symbols denote the measurement 
of energy fluctuations. 



show no divergent behaviour. This is consistent with a KT-like transition and makes a 
first- or second-order transition unlikely. 

The magnetic susceptibility was calculated from fluctuations in the magnetization per 
point 

x = L 2 (m 2 ), rh = ^-^s(x,y) (11) 

x,y 

in the conventional manner and using improved estimators. Finite-size effects were studied 
by comparing data on different lattice sizes. We find that for correlation lengths £ < L/7 
no lattice size dependence of x within the statistical errors can be seen. 

In the high-temperature phase the correlation length was extracted from the 'zero 
momentum' spin-spin correlation function 

r,(x) = (s(x) • 5(o)) , s(x) = yE s>> v) (12) 

L y 

by fitting the data with a single cosh or cosh + const in the interval Xq < x < L/2. Most 
of the time, constant contributions were negligible, because they are of order exp(—L/^). 
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To determine the influence of 'excitations' (smaller eigenvalues of the transfer matrix), 
we compared the results for different minimal distance xq. The correlations are always 
dominated by the lowest state of the transfer operator 'Hamiltonian', so that we had 
to omit only few points. As for the susceptibility we used conventional and improved 
estimators. The latter have substantially smaller statistical errors, because only positive 
terms contribute to the correlation function. For (3 > /3 C it is expected that the long 
distance behaviour of the spin-spin correlation function is given by a power-law decrease 
with a critical exponent rj which is a function of ft, but we made no investigations to this 
point. 



2.2 Numerical results in the high-temperature phase 

In the following we will analyze the critical behaviour of £ and \. Therefore, we have 
performed least square fits with three different forms (according to eqs. (|3|) and ([?])). In 
the first case we used a four-parameter Kosterlitz-Thouless fit (KT4), i.e. 

ln(f) = a + bt~ u , (13) 

with (3 C = f3 + 1 as the fourth parameter. In the second case the value of v was fixed to 0.5 
(KT3). Assuming a power- law behaviour of the divergence we obtained the third case: 

ln(£) = a - v ln(t) . (14) 

Systematic errors had been estimated by varying the range of points to be fitted and by 
replacing t — j3 c — (3 by t = T — T c . We used data subsets consisting of all points with 
6 < f < L/7 and 12 < f < L/7, respectively!]. This corresponds to 0.67 < f3 < 0.78 (20 
points) and 0.72 < (3 < 0.78 (15 points). Errors for the fitted parameters were computed 
by performing fits on data sets to a Gaussian distribution of £(/3) with variance A£(/3). 
All calculations were performed correspondingly for \. 

Our results are summarized in table [l|. Additionally one example is plotted in fig. [2|. 
One observes that the KT like fits look consistent (except one of the KT3 fits for £(/?)) 
and that both types always give a similar \ 2 per degree of freedom (dof). The influence 
of the lower bound of the fitted interval is small. This situation changes for the power-law 
fits. The different intervals yield large changes in the fitted parameter. Also, PL fits 
have always larger x 2 /dof than the KT fits. The data for the correlation length make a 
power-law divergence very unlikely. 

Taking v as a free parameter (KT4) results in large errors. This situation is similar as 
in the XY model fl~3| , [14]. The reason is that in the four-parameter space a valley for x 2 
exists, which is narrow in three dimensions, but flat in the fourth. The different fits with 
v as a free parameter yield the estimate 

v = 0.46(7) . (15) 

Holding v or T c fixed leads to a clear minimum and much smaller errors in the remaining 
parameters. From the /3 c 's of the three-parameter Kosterlitz-Thouless fits and the power- 
law fits we got the following values: 

^KT3 = 0.8354(20) , (16) 

1 There is one exception: for the point (3 = 0.7725 we also got L/£ « 480/69.2 w 6.94. 
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Figure 2: The curves shown are the best fit for a Kosterlitz-Thouless be- 
haviour with fixed v (dotted) and a power-law ansatz (dashed). We used all 
correlation lengths £ > 6. Statistical errors are too small for a visualization. 
The critical values of /3 are visualized by vertical lines. 



&, PL = 0.8079(30) . (17) 

We further analyzed our data by investigating the relation between \ and £. To check 
the validity of eqn. (|5|) we plot ln(x/£' ) versus Zn(£). For the predicted KT value 77 = 1/4 
we should see a horizontal line. A different value of r] would correspond to a straight line 
of non-zero slope. Indeed, fig. |3] shows a negative slope, but with decreasing absolute 
value for increasing £. Performing straight line fits for the points £ > 12 results in an 
estimate of 

77 = 0.270(3) . (18) 

For the points £ > 50 we got 77 = 0.259(17). The difference to the theoretical value can 
perhaps be explained by logarithmical corrections as discussed later. 

It should be noted that the values of the critical exponent v = 0.48(10) and 77 = 0.267 
of the XY model in Villain's formulation |14| (calculated with same methods) coincide 
with our grain boundary model within statistical errors. Also, the behaviour of /n(x/£ 7 / 4 ) 
as a function of Zn(£ ) (fig. 0) is qualitatively the same in both cases. Hence it looks likely 
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Table 1: Best fit parameters for the critical behaviour of the correla- 
tion length and susceptibility. For £ > 6 we have fitted 20 points in 
0.67 < (3 < 0.78. In the case £ > 12 we used 15 points in the range 
0.72 < (3 < 0.78. 







KT4 


KT3 


PL 






£ > 6 


£ > 12 


£ > 6 


£ > 12 


£ > 6 


£ > 12 




m(aj 


-3.62(33) 


-2.92(88) 


-2.08(1) 


-2.01(2) 


0.176(1) 


0.148(2) 




b 


2.88(28) 


2.30(73) 


1.62(1) 


1.59(1) 








I3 C 


0.8295(12) 


0.8323(36) 


0.8382(2) 


0.8373(3) 


0.8097(1) 


0.8118(2) 




V 


0.350(33) 


0.403(71) 


0.5 


0.5 


1.816(3) 


1.901(6) 




X /dot 


1.10 


1.17 


3.37 


1.22 


28.41 


6.59 




m(a) 


-1.45(13) 


-1.25(40) 


-1.32(1) 


-1.32(2) 


0.713(2) 


0.594(4) 




b 


1.84(11) 


1.67(33) 


1.73(1) 


1.73(1) 
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is 


0.480(19) 


0.510(62) 


0.5 


0.5 


1.557(3) 


1.669(5) 




X 2 /dof 


1.00 


1.16 


1.02 


1.07 


80.64 


14.08 




ln(a) 


-3.7(12) 


-1.0(20) 


-2.46(4) 


-2.32(10) 


0.187(4) 


0.126(9) 




b 


3.7(10) 


1.7(15) 


2.75(2) 


2.69(5) 






x(P) 


13c 


0.8318(43) 


0.8335(93) 


0.8372(4) 


0.8361(9) 


0.8069(3) 


0.8099(6) 




1/(7) 


0.42(7) 


0.64(17) 


0.5 


0.5 


2.974(12) 


3.169(34) 




X 2 /dof 


0.94 


1.25 


0.97 


1.27 


4.84 


2.49 




ln(a) 


-0.8(5) 


0.5(15) 


-1.12(3) 


-1.13(8) 


1.98(2) 


1.33(5) 




b 


2.6(4) 


1.6(12) 


2.91(2) 


2.91(4) 






X(T) 


(3c 


0.8359(40) 


0.8455(99) 


0.8332(4) 


0.8333(8) 


0.8024(3) 


0.8066(6) 




1/(7) 


0.54(6) 


0.71(27) 


0.5 


0.5 


2.524(9) 


2.771(26) 




x 2 /dof 


0.93 


1.25 


0.91 


1.23 


11.88 


3.06 



that the grain boundary model is in the same universality class as the XY model. Also, 
a detailed comparison indicates a faster convergence to the KT behaviour for the lattice 
grain boundary model than for the XY model in Villain's formulation. The same analysis 
for the XY model in the cosine form yields non-unique results, since the data of |13|, [15| 
have larger statistical errors. 

2.3 Numerical results at the transition point 

We now come to the simulations near the transition point. In the following we will 
concentrate on the fourth-order cumulant 



3{m z ) z 



and the susceptibility. At the transition point finite-size scaling (FSS) implies scale in- 
variance of U and x ~ L 2 ~ v for large enough lattices. We use this FSS behaviour to locate 
f3 c . For these simulations we used additional lattices of size L 2 = 36 2 . 



7 



>< 



1.15 
1.1 
1.05 
1 

0.95 
0.9 



i i i i i i i i i i i 



60 



120 



240 



i | i 
480 



i i | 
960 



i i i i i i i i i i i 



5 % 



i i i : 



8 



i i i 



i i i i 



1 2 3 4 5 
ln(f) 

Figure 3: Test of the scaling relation \ ~ £ 7 ^ 4 - The tendency to decrease 
implies r\ > 1/4. 



At the transition point the susceptibility should diverge with the size of the system as 

X ~ L 2 ^ . (20) 

For (3 > (3 C 7] is a function of the temperature. For /3's below the transition point, one 
has to take corrections for finite correlation lengths of the order 0(L/£) into account. In 
fig. |3] we plot ln(x/ ' L 7 / 4 ) versus ln(L) for three different /3's near the transition point. The 
slope gives the deviation from rj = 1/4. We extracted the values of rj from linear fits of 
the asymptotic behaviour. Our results are collected in table 0. Obviously, a requirement 
of rj = 0.25 yields (3 C ~ 0.826. Since there is a tendency of the slope to decrease with 



Table 2: The exponent r] for different /3's as obtained from FSS. 



p 


j] 


X 2 /dof 


used lattices 


0.82 


0.2663(13) 


0.89 


L — 120 - 480 


0.825 


0.2535(10) 


2.91 


L — 120 - 960 


0.83 


0.2404(7) 


0.01 


L = 60 - 480 


0.84 


0.2242(13) 


1.49 


L = 60 - 240 
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Figure 4: Finite-size scaling of the susceptibility for various inverse tem- 
peratures. The slope gives the deviation from rj — 1/4. 



larger lattices, the real scale invariance may take place at larger inverse temperatures. On 
the other hand, with /3 c ,kt3 ~ 0.836 we got an 77 which is about 8 % below the theoretical 



value. Again this situation is similar as in the XY model [13], [14 . 

Our results for the dependence of U on the lattice size are shown in fig. |. Apart from 
the fact that the statistical errors are larger, we get a similar estimate as before. Also the 
slopes of Zn(x/L 7 / 4 ) and of U both decrease with increasing L. Therefore, the value of (3 C 
can be taken as a lower bound of the large L limit. If we compare this value with those 
of the KT and PL fits (eqn. ( JIB] ) and QT7D) we find that /3 c ,kt is in better agreement than 
/3c,pl- This is an additional hint favouring a KT transition. 

2.4 Logarithmic corrections to KT scaling laws 

Assuming an exponential divergence as predicted by Kosterlitz and Thouless, one can ask 
for corrections of the scaling behaviour of £ and x- The renormalization group analysis 
of KT yields the following multiplicative corrections to eqn. (|5|): 

X ~ t r e~\ r = ~. (21) 

lb 

In the following we will analyze the data under the assumption that v — 1/2, r] — 1/4 
and determine if logarithmical corrections can explain the appearing discrepancies from 
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Figure 5: Finite-size scaling of the cumulant in the vicinity of the transition 
point. 



the leading scaling behaviour. Therefore in fig. ^ we plotted ln{x/£?^) as a function of 
—ln(t) = —ln(/3 c — P), where we have taken f3 c = 0.8354. The data are consistent with a 
straight line assumption, but from the slope (for the points £ > 6) we obtained 



0.070(5) 



(22) 



which is completely different from the theoretical value. For the interval £ > 12 we got 
r = 0.056(9). This means that the value of r decreases, so that we probably have not 
reached the scaling region. On the other hand a different f3 c results in only small changes 
of r. 

Alternatively one can also try to replace —ln(t) by 2ln(ln(£)) in fig. ^|, as it follows 
from eqn. (|3]). This has the advantage, that no information about (3 C is needed. It is 
clear that the relation is only fulfilled in the limit t — > 0, where h^t~ u ^> \ln(a^)\. For 
the actual values of (3 this is not the case, as can be calculated using the parameters of 
table |l| Nevertheless one gets a linear behaviour which leads to r = 0.043(3) (£ > 6) and 
r = 0.036(6) (£ > 12), respectively. All data are summarized in table || 

A different approach to estimate r is based one FSS. At the transition point (£ 3> L) 
one expects 
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(ln(L)) 
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Figure 6: Logarithmic corrections to the scaling of \ ~ C^ A - The slope 
yields — r. The dashed line results from a fit of all points with £ > 6, i.e. 
-Zn(/3 C - 0) > 1.77. 



In fig. |7| we show the data at (3 — 0.83 ~ (3 C . Indeed we can observe the predicted linear 
behaviour. From the slope we obtain 

r = -0.0233(10) . (24) 

Although the value is closer to the theoretical one, it is still in disagreement with it. A 
higher (lower) value of (3 would result in a decrease (increase) of r. At /3 = 0.84 we got 
r = -0.059(4). 

The analysis for multiplicative corrections in the XY model yields similar results |[T6|, 



I7||, i.e. positive values of the order (9(0.05) from the data of the high-temperature phase 



and a negative value of approximate —0.03 from finite-size scaling. 



3 Conclusions 

We presented a detailed Monte Carlo study of a two-dimensional grain boundary model 
on the lattice. The investigations were performed in the high-temperature phase and near 
the phase transition point. 

The behaviour of the specific heat as a function of the inverse temperature f3 was char- 
acterized by a smooth peak at an inverse temperature below the critical value. Moreover 
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Table 3: Values of the logarithmic correction exponent r obtained from 
different methods. 



method 


r 


X Vdof 


points 


H X /e /4 ) vs. ln(/3 c - p) 


f > 6 


0.070(5) 


0.86 


20 


£> 12 


0.056(9) 


0.57 


15 


Hx/e /4 ) vs. ln(ln(0) 


£ > 6 


0.043(3) 


0.67 


20 


£> 12 


0.036(6) 


0.56 


15 


ln(x/^ 7/4 ) vs. ln(ln(L)) 


/3 = 0.83 


-0.0233(10) 


0.26 


5 


P = 0.84 


-0.059(4) 


4.20 


3 



for lattices of L = 120 or larger no finite-size effects were seen. This is consistent with 
the KT scenario. 

In the high-temperature phase we examined the dependency of the correlation length 
and susceptibility on f3 and T. We showed that the data are always in agreement with a 
KT-like divergence. The critical exponent v was estimated as v = 0.46(7). In all cases fits 
with a power-law divergence are characterized by a larger % 2 /dof than the corresponding 
KT fit. The exponent i] was derived from the relation of £ and x for /3 — > f3 c . We got 
7} = 0.270(3). 

The simulations in the vicinity of the transition point were used to measure the finite- 
size scaling of the fourth-order cumulant and susceptibility. The resulting estimates of (3 C 
are closer to the value from the KT fits than to the estimate from the PL fits. Deviations 
may be explained by the fact that the lattices used are still too small. 

Additionally we discussed an attempt to explain differences of the critical exponent rj 
from the predicted KT values by multiplicative corrections. We showed that both finite- 
size scaling at f3 c and simulations in the high-temperature phase are consistent with such 
an ansatz. However the values of r are inconsistent and far from the theoretical value 
r = —1/16. As before the reason might be too small lattice sizes. 

A comparison of our results showed that the data are in accordance with the ones 
of the XY model. The universal critical exponents agree within statistical errors, while 
multiplicative corrections are qualitatively the same. We take this as evidence that both 
models lie in the same universality class. 
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